
# ------ Libraries ------ #
library(texreg)
library(xergm)
library(data.table)
library(ggplot2)
set.seed(1912)

setwd("YOUR WORKING DIRECTORY")


# --- Notification:

note <- function(message=message) {system(sprintf("/usr/local/bin/terminal-notifier -title 'RStudio' -message '%s' -activate org.rstudio.RStudio",
                                                  message), ignore.stdout=TRUE, ignore.stderr=TRUE, wait=FALSE)}







# --- In sample GOF --- #

mid.nw_full <- readRDS("data/pauls_cranmer/prepped/system/mid.nw_full.rds")
directcont.nw_full <- readRDS("data/pauls_cranmer/prepped/system/directcont.nw_full.rds")
capratio_full <- readRDS("data/pauls_cranmer/prepped/system/capratio_full.rds")
tradedep_full <- readRDS("data/pauls_cranmer/prepped/system/tradedep_full.rds")
igosec_full <- readRDS("data/pauls_cranmer/prepped/system/igosec_full.rds")
igoecon_full <- readRDS("data/pauls_cranmer/prepped/system/igoecon_full.rds")





# --- Model 1:
mod.1 <- readRDS("inference/system/model.1_baseline.rds")
gof1 <- gof(mod.1, nsim = 100);note("Done.")
#saveRDS(gof1, file = "inference/gof/gof1.rds")



# --- Model 2:
weak_comms <- readRDS("data/comm_detection/detected_comms/weak/weak_ties_50d_15x.rds")[1970:1989] 
mod.2 <- readRDS("inference/system/weak/model.2_weak_weak_ties_50d_15x.rds")
gof2 <- gof(mod.2, nsim = 100);note("Done.")
#saveRDS(gof2, file = "inference/gof/gof2.rds")



# --- Model 3:
strong_comms <- readRDS("data/comm_detection/detected_comms/strong/strong_ties_masternodes_prop0.2.rds")[1971:1990] 
mod.3 <- readRDS("inference/system/strong/model.3_strong_strong_ties_masternodes_prop0.2.rds")
gof3 <- gof(mod.3, nsim = 100);note("Done.")
#saveRDS(gof3, file = "inference/gof/gof3.rds")



# --- Model 4:
strong_comms <- readRDS("data/comm_detection/detected_comms/strong/strong_ties_masternodes_prop0.2.rds")[1971:1990] 
mod.4 <- readRDS("inference/system/strong/model.4_strong_nocontig_strong_ties_masternodes_prop0.2.rds")
gof4 <- gof(mod.4, nsim = 100);note("Done.")
#saveRDS(gof4, file = "inference/gof/gof4.rds")



# --- Model 5:
directcont_sub_w <- readRDS("data/pauls_cranmer/prepped/subsets/weak_nodes/predictors/directcont_sub_w.rds")
capratio_sub_w <- readRDS("data/pauls_cranmer/prepped/subsets/weak_nodes/predictors/capratio_sub_w.rds")
tradedep_sub_w <- readRDS("data/pauls_cranmer/prepped/subsets/weak_nodes/predictors/tradedep_sub_w.rds")
igosec_sub_w <- readRDS("data/pauls_cranmer/prepped/subsets/weak_nodes/predictors/igosec_sub_w.rds")
igoecon_sub_w <- readRDS("data/pauls_cranmer/prepped/subsets/weak_nodes/predictors/igoecon_sub_w.rds")
mid_sub_w <- readRDS("data/pauls_cranmer/prepped/subsets/weak_nodes/outcome_nets/mid_sub_w_50d_15x.rds")

mod.5 <- readRDS("inference/community_subset/weak/model.5_mid_sub_weak_mid_sub_w_50d_15x.rds")
gof5 <- gof(mod.5, nsim = 100);note("Done.")
#saveRDS(gof5, file = "inference/gof/gof5.rds")



# --- Model  6:
directcont_sub_s <- readRDS("data/pauls_cranmer/prepped/subsets/strong_nodes/predictors/nodeset_weak/directcont_sub_s_weaknodes.rds")
capratio_sub_s <- readRDS("data/pauls_cranmer/prepped/subsets/strong_nodes/predictors/nodeset_weak/capratio_sub_s_weaknodes.rds")
tradedep_sub_s <- readRDS("data/pauls_cranmer/prepped/subsets/strong_nodes/predictors/nodeset_weak/tradedep_sub_s_weaknodes.rds")
igosec_sub_s <- readRDS("data/pauls_cranmer/prepped/subsets/strong_nodes/predictors/nodeset_weak/igosec_sub_s_weaknodes.rds")
igoecon_sub_s <- readRDS("data/pauls_cranmer/prepped/subsets/strong_nodes/predictors/nodeset_weak/igoecon_sub_s_weaknodes.rds")
mid_sub_s <- readRDS("data/pauls_cranmer/prepped/subsets/strong_nodes/outcome_nets/mid_sub_s_weaknodes_prop0.2.rds")

mod.6 <- readRDS("inference/community_subset/strong/model.6_mid_sub_strong_mid_sub_s_weaknodes_prop0.2.rds")
gof6 <- gof(mod.6, nsim = 100);note("Done.")
#saveRDS(gof6, file = "inference/gof/gof6.rds")





# --- GOF Plots:
gof1 <- readRDS("inference/gof/gof1.rds")
gof2 <- readRDS("inference/gof/gof2.rds")
gof3 <- readRDS("inference/gof/gof3.rds")
gof4 <- readRDS("inference/gof/gof4.rds")
gof5 <- readRDS("inference/gof/gof5.rds")
gof6 <- readRDS("inference/gof/gof6.rds")

pdf("inference/gof/m1_deg_gof.pdf", width = 12, height = 8)
par(mar=c(5.1,4.1,1,2.1));par(bty="n",lwd=.2, cex=2, lty=1)
btergm::plot(gof1$Degree, boxplot.lwd=1.8, transform = function(x) log(x + 1), main="");dev.off()

pdf("inference/gof/m1_esp_gof.pdf", width = 12, height = 8)
par(mar=c(5.1,4.1,1,2.1));par(bty="n",lwd=.2, cex=2, lty=1)
btergm::plot(gof1$`Edge-wise shared partners`, boxplot.lwd=1.8, transform = function(x) log(x + 1), main="Model 1");dev.off()

pdf("inference/gof/m1_mod_gof.pdf", width = 12, height = 8)
par(mar=c(5.1,4.1,1,2.1));par(bty="n",lwd=.2, cex=2, lty=1)
btergm::plot(gof1$`Modularity (walktrap)`, boxplot.lwd=1.8, main = "");dev.off()

pdf("inference/gof/m2_deg_gof.pdf", width = 12, height = 8)
par(mar=c(5.1,4.1,1,2.1));par(bty="n",lwd=.2, cex=2, lty=1)
btergm::plot(gof2$Degree, boxplot.lwd=1.8, transform = function(x) log(x + 1), main="");dev.off()

pdf("inference/gof/m2_esp_gof.pdf", width = 12, height = 8)
par(mar=c(5.1,4.1,1,2.1));par(bty="n",lwd=.2, cex=2, lty=1)
btergm::plot(gof2$`Edge-wise shared partners`, boxplot.lwd=1.8, transform = function(x) log(x + 1), main="Model 2");dev.off()

pdf("inference/gof/m2_mod_gof.pdf", width = 12, height = 8)
par(mar=c(5.1,4.1,1,2.1));par(bty="n",lwd=.2, cex=2, lty=1)
btergm::plot(gof2$`Modularity (walktrap)`, boxplot.lwd=1.8, main = "");dev.off()

pdf("inference/gof/m3_deg_gof.pdf", width = 12, height = 8)
par(mar=c(5.1,4.1,1,2.1));par(bty="n",lwd=.2, cex=2, lty=1)
btergm::plot(gof3$Degree, boxplot.lwd=1.8, transform = function(x) log(x + 1), main="");dev.off()

pdf("inference/gof/m3_esp_gof.pdf", width = 12, height = 8)
par(mar=c(5.1,4.1,1,2.1));par(bty="n",lwd=.2, cex=2, lty=1)
btergm::plot(gof3$`Edge-wise shared partners`, boxplot.lwd=1.8, transform = function(x) log(x + 1), main="Model 3");dev.off()

pdf("inference/gof/m3_mod_gof.pdf", width = 12, height = 8)
par(mar=c(5.1,4.1,1,2.1));par(bty="n",lwd=.2, cex=2, lty=1)
btergm::plot(gof3$`Modularity (walktrap)`, boxplot.lwd=1.8, main = "");dev.off()

pdf("inference/gof/m4_deg_gof.pdf", width = 12, height = 8)
par(mar=c(5.1,4.1,1,2.1));par(bty="n",lwd=.2, cex=2, lty=1)
btergm::plot(gof4$Degree, boxplot.lwd=1.8, transform = function(x) log(x + 1), main="");dev.off()

pdf("inference/gof/m4_esp_gof.pdf", width = 12, height = 8)
par(mar=c(5.1,4.1,1,2.1));par(bty="n",lwd=.2, cex=2, lty=1)
btergm::plot(gof4$`Edge-wise shared partners`, boxplot.lwd=1.8, transform = function(x) log(x + 1), main="Model 4");dev.off()

pdf("inference/gof/m4_mod_gof.pdf", width = 12, height = 8)
par(mar=c(5.1,4.1,1,2.1));par(bty="n",lwd=.2, cex=2, lty=1)
btergm::plot(gof4$`Modularity (walktrap)`, boxplot.lwd=1.8, main = "");dev.off()

pdf("inference/gof/m5_deg_gof.pdf", width = 12, height = 8)
par(mar=c(5.1,4.1,1,2.1));par(bty="n",lwd=.2, cex=2, lty=1)
btergm::plot(gof5$Degree, boxplot.lwd=1.8, transform = function(x) log(x + 1), main="");dev.off()

pdf("inference/gof/m5_esp_gof.pdf", width = 12, height = 8)
par(mar=c(5.1,4.1,1,2.1));par(bty="n",lwd=.2, cex=2, lty=1)
btergm::plot(gof5$`Edge-wise shared partners`, boxplot.lwd=1.8, transform = function(x) log(x + 1), main="Model 5");dev.off()

pdf("inference/gof/m5_mod_gof.pdf", width = 12, height = 8)
par(mar=c(5.1,4.1,1,2.1));par(bty="n",lwd=.2, cex=2, lty=1)
btergm::plot(gof5$`Modularity (walktrap)`, boxplot.lwd=1.8, main = "");dev.off()

pdf("inference/gof/m6_deg_gof.pdf", width = 12, height = 8)
par(mar=c(5.1,4.1,1,2.1));par(bty="n",lwd=.2, cex=2, lty=1)
btergm::plot(gof6$Degree, boxplot.lwd=1.8, transform = function(x) log(x + 1), main="");dev.off()

pdf("inference/gof/m6_esp_gof.pdf", width = 12, height = 8)
par(mar=c(5.1,4.1,1,2.1));par(bty="n",lwd=.2, cex=2, lty=1)
btergm::plot(gof6$`Edge-wise shared partners`, boxplot.lwd=1.8, transform = function(x) log(x + 1), main="Model 6");dev.off()

pdf("inference/gof/m6_mod_gof.pdf", width = 12, height = 8)
par(mar=c(5.1,4.1,1,2.1));par(bty="n",lwd=.2, cex=2, lty=1)
btergm::plot(gof6$`Modularity (walktrap)`, boxplot.lwd=1.8, main = "");dev.off()













# --- Out of sample predictions --- #


# --- Model 1 (Baseline):

gof.list.1 <- list()
for (i in 1:15){

  print(i)
  
  j=i+4
  k=i+5
  lst=i:j
  
  
  directcont.nw.gof1 <- directcont.nw_full[lst]
  capratio.gof1 <- capratio_full[lst]
  tradedep.gof1 <- tradedep_full[lst]
  igosec.gof1 <- igosec_full[lst]
  igoecon.gof1 <- igoecon_full[lst]
  
  set.seed(1912)
  model.1.Mult.gof <- btergm(
    mid.nw_full[lst] ~ 
      edges + 
      altkstar(2, fixed = TRUE) + 
      cycle(4) + 
      gwesp(0, fixed = TRUE) +
      nodemix("democ", -3) +
      edgecov(directcont.nw.gof1) +
      edgecov(capratio.gof1) + 
      edgecov(tradedep.gof1) + 
      edgecov(igosec.gof1) + 
      edgecov(igoecon.gof1), 
    R = 1000
    , parallel = "multicore", ncpus = 3
  )
  
  
  directcont.nw.gof1 <- directcont.nw_full[k]
  capratio.gof1 <- capratio_full[k]
  tradedep.gof1 <- tradedep_full[k]
  igosec.gof1 <- igosec_full[k]
  igoecon.gof1 <- igoecon_full[k]
  
  set.seed(1912)
  gof.ofs.E1 <- gof(model.1.Mult.gof, nsim=50, target = mid.nw_full[[k]],
                    formula = mid.nw_full[k] ~ 
                      edges + 
                      altkstar(2, fixed = TRUE) + 
                      cycle(4) + 
                      gwesp(0, fixed = TRUE) + 
                      nodemix("democ", -3) + 
                      edgecov(directcont.nw.gof1) + 
                      edgecov(capratio.gof1) + 
                      edgecov(tradedep.gof1) + 
                      edgecov(igosec.gof1) + 
                      edgecov(igoecon.gof1),
                    coef = coef(model.1.Mult.gof), statistics = c(rocpr), 
                    parallel = "multicore", ncpus = 3)
  
  gof.list.1[[i]] <- gof.ofs.E1
  
};note("Done.")

#saveRDS(gof.list.1, file="inference/gof/oos_model.1.rds")










# --- Model 2 (Baseline + Weak Comms)
weak_comms <- readRDS("data/comm_detection/detected_comms/weak/weak_ties_50d_15x.rds")[1970:1989] 
weak_comms <- adjust(weak_comms, mid.nw_full, remove = TRUE, add = TRUE, value = 0)


gof.list.2 <- list()
for (i in 1:15){
  
  print(i)
  
  j=i+4
  k=i+5
  lst=i:j
  
  
  directcont.nw.gof1 <- directcont.nw_full[lst]
  capratio.gof1 <- capratio_full[lst]
  tradedep.gof1 <- tradedep_full[lst]
  igosec.gof1 <- igosec_full[lst]
  igoecon.gof1 <- igoecon_full[lst]
  weak_comms.gof1 <- weak_comms[lst]
  
  set.seed(1912)
  model.1.Mult.gof <- btergm(
    mid.nw_full[lst] ~ 
      edges + 
      altkstar(2, fixed = TRUE) + 
      cycle(4) + 
      gwesp(0, fixed = TRUE) + 
      nodemix("democ", -3) + 
      edgecov(directcont.nw.gof1) + 
      edgecov(capratio.gof1) + 
      edgecov(tradedep.gof1) + 
      edgecov(igosec.gof1) + 
      edgecov(igoecon.gof1) +
      edgecov(weak_comms.gof1), 
    R = 1000
    , parallel = "multicore", ncpus = 3
  )
  
  
  directcont.nw.gof1 <- directcont.nw_full[k]
  capratio.gof1 <- capratio_full[k]
  tradedep.gof1 <- tradedep_full[k]
  igosec.gof1 <- igosec_full[k]
  igoecon.gof1 <- igoecon_full[k]
  weak_comms.gof1 <- weak_comms[k]
  
  
  set.seed(1912)
  gof.ofs.E1 <- gof(model.1.Mult.gof, nsim=50, target = mid.nw_full[[k]],
                    formula = mid.nw_full[k] ~ 
                      edges + 
                      altkstar(2, fixed = TRUE) + 
                      cycle(4) + 
                      gwesp(0, fixed = TRUE) + 
                      nodemix("democ", -3) + 
                      edgecov(directcont.nw.gof1) + 
                      edgecov(capratio.gof1) + 
                      edgecov(tradedep.gof1) + 
                      edgecov(igosec.gof1) + 
                      edgecov(igoecon.gof1) +
                      edgecov(weak_comms.gof1),
                    coef = coef(model.1.Mult.gof), statistics = c(rocpr), 
                    parallel = "multicore", ncpus = 3)
  
  gof.list.2[[i]] <- gof.ofs.E1
  
};note("Done.")

  
#saveRDS(gof.list.2, file="inference/gof/oos_model.2.rds")











# --- Model 3 (Baseline + Strong Comms)
strong_comms <- readRDS("data/comm_detection/detected_comms/strong/strong_ties_masternodes_prop0.2.rds")[1971:1990]
strong_comms <- adjust(strong_comms, mid.nw_full, remove = TRUE, add = TRUE, value = 0)


gof.list.3 <- list()
for (i in 1:15){
  
  print(i)
  
  j=i+4
  k=i+5
  lst=i:j
  
  
  directcont.nw.gof1 <- directcont.nw_full[lst]
  capratio.gof1 <- capratio_full[lst]
  tradedep.gof1 <- tradedep_full[lst]
  igosec.gof1 <- igosec_full[lst]
  igoecon.gof1 <- igoecon_full[lst]
  strong_comms.gof1 <- strong_comms[lst]
  
  set.seed(1912)
  model.1.Mult.gof <- btergm(
    mid.nw_full[lst] ~ 
      edges + 
      altkstar(2, fixed = TRUE) + 
      cycle(4) + 
      gwesp(0, fixed = TRUE) + 
      nodemix("democ", -3) + 
      edgecov(directcont.nw.gof1) + 
      edgecov(capratio.gof1) + 
      edgecov(tradedep.gof1) + 
      edgecov(igosec.gof1) + 
      edgecov(igoecon.gof1) +
      edgecov(strong_comms.gof1), 
    R = 1000
    , parallel = "multicore", ncpus = 3
  )
  
  
  directcont.nw.gof1 <- directcont.nw_full[k]
  capratio.gof1 <- capratio_full[k]
  tradedep.gof1 <- tradedep_full[k]
  igosec.gof1 <- igosec_full[k]
  igoecon.gof1 <- igoecon_full[k]
  strong_comms.gof1 <- strong_comms[k]
  
  set.seed(1912)
  gof.ofs.E1 <- gof(model.1.Mult.gof, nsim=50, target = mid.nw_full[[k]],
                    formula = mid.nw_full[k] ~ 
                      edges + 
                      altkstar(2, fixed = TRUE) + 
                      cycle(4) + 
                      gwesp(0, fixed = TRUE) + 
                      nodemix("democ", -3) + 
                      edgecov(directcont.nw.gof1) + 
                      edgecov(capratio.gof1) + 
                      edgecov(tradedep.gof1) + 
                      edgecov(igosec.gof1) + 
                      edgecov(igoecon.gof1) +
                      edgecov(strong_comms.gof1),
                    coef = coef(model.1.Mult.gof), statistics = c(rocpr), 
                    parallel = "multicore", ncpus = 3)
  
  gof.list.3[[i]] <- gof.ofs.E1
  
};note("Done.")


#saveRDS(gof.list.3, file="inference/gof/oos_model.3.rds")











# --- Model 4 (Baseline + Strong Comms -- no contig)

gof.list.4 <- list()
for (i in 1:15){
  
  print(i)
  
  j=i+4
  k=i+5
  lst=i:j
  
  capratio.gof1 <- capratio_full[lst]
  tradedep.gof1 <- tradedep_full[lst]
  igosec.gof1 <- igosec_full[lst]
  igoecon.gof1 <- igoecon_full[lst]
  strong_comms.gof1 <- strong_comms[lst]
  
  set.seed(1912)
  model.1.Mult.gof <- btergm(
    mid.nw_full[lst] ~ 
      edges + 
      altkstar(2, fixed = TRUE) + 
      cycle(4) + 
      gwesp(0, fixed = TRUE) + 
      nodemix("democ", -3) + 
      edgecov(capratio.gof1) + 
      edgecov(tradedep.gof1) + 
      edgecov(igosec.gof1) + 
      edgecov(igoecon.gof1) +
      edgecov(strong_comms.gof1), 
    R = 1000
    , parallel = "multicore", ncpus = 3
  )
  
  
  capratio.gof1 <- capratio_full[k]
  tradedep.gof1 <- tradedep_full[k]
  igosec.gof1 <- igosec_full[k]
  igoecon.gof1 <- igoecon_full[k]
  strong_comms.gof1 <- strong_comms[k]
  
  set.seed(1912)
  gof.ofs.E1 <- gof(model.1.Mult.gof, nsim=50, target = mid.nw_full[[k]],
                    formula = mid.nw_full[k] ~ 
                      edges + 
                      altkstar(2, fixed = TRUE) + 
                      cycle(4) + 
                      gwesp(0, fixed = TRUE) + 
                      nodemix("democ", -3) + 
                      edgecov(capratio.gof1) + 
                      edgecov(tradedep.gof1) + 
                      edgecov(igosec.gof1) + 
                      edgecov(igoecon.gof1) +
                      edgecov(strong_comms.gof1),
                    coef = coef(model.1.Mult.gof), statistics = c(rocpr), 
                    parallel = "multicore", ncpus = 3)
  
  gof.list.4[[i]] <- gof.ofs.E1
  
};note("Done.")

#saveRDS(gof.list.4, file="inference/gof/oos_model.4.rds")











# --- Model 5 (Network effects + node effects for weak comms)
# NB: because bridging ties are not present in every year, we won't be able to estimate a coefficient during those absent years.
# so, either fit and estimate out-of-sample accuracy without the nodefactor("weak_bridge") term or subset the lists to
# only those years in which weak bridges are present. the results are substantively similar either way, so the models below 
# opt to retain the entire time range and skip the nodefactor("weak_bridge") term.

mid_sub_w <- readRDS("data/pauls_cranmer/prepped/subsets/weak_nodes/outcome_nets/mid_sub_w_100d_15x.rds")
directcont_sub_w <- readRDS("data/pauls_cranmer/prepped/subsets/weak_nodes/predictors/directcont_sub_w.rds")
capratio_sub_w <- readRDS("data/pauls_cranmer/prepped/subsets/weak_nodes/predictors/capratio_sub_w.rds")
tradedep_sub_w <- readRDS("data/pauls_cranmer/prepped/subsets/weak_nodes/predictors/tradedep_sub_w.rds")
igosec_sub_w <- readRDS("data/pauls_cranmer/prepped/subsets/weak_nodes/predictors/igosec_sub_w.rds")
igoecon_sub_w <- readRDS("data/pauls_cranmer/prepped/subsets/weak_nodes/predictors/igoecon_sub_w.rds")

gof.list.5 <- list()
for (i in 1:15){

  print(i)
  
  j=i+4
  k=i+5
  lst=i:j
  
  
  directcont.nw.gof1 <- directcont_sub_w[lst]
  capratio.gof1 <- capratio_sub_w[lst]
  tradedep.gof1 <- tradedep_sub_w[lst]
  igosec.gof1 <- igosec_sub_w[lst]
  igoecon.gof1 <- igoecon_sub_w[lst]

  set.seed(1912)
  model.1.Mult.gof <- btergm(
    mid_sub_w[lst] ~ 
      edges + 
      altkstar(2, fixed = TRUE) + 
      cycle(4) + 
      gwesp(0, fixed = TRUE) + 
      edgecov(directcont.nw.gof1) +
      edgecov(capratio.gof1) + 
      edgecov(tradedep.gof1) + 
      edgecov(igosec.gof1) + 
      edgecov(igoecon.gof1) +
      nodematch("weak_comm"), 
    R = 1000
    , parallel = "multicore", ncpus = 3
  )
  
  
  directcont.nw.gof1 <- directcont_sub_w[k]
  capratio.gof1 <- capratio_sub_w[k]
  tradedep.gof1 <- tradedep_sub_w[k]
  igosec.gof1 <- igosec_sub_w[k]
  igoecon.gof1 <- igoecon_sub_w[k]
  
  set.seed(1912)
  gof.ofs.E1 <- gof(model.1.Mult.gof, nsim=50, target = mid_sub_w[[k]],
                    formula = mid_sub_w[k] ~ 
                      edges + 
                      altkstar(2, fixed = TRUE) + 
                      cycle(4) + 
                      gwesp(0, fixed = TRUE) + 
                      edgecov(directcont.nw.gof1) +
                      edgecov(capratio.gof1) + 
                      edgecov(tradedep.gof1) + 
                      edgecov(igosec.gof1) + 
                      edgecov(igoecon.gof1) +
                      nodematch("weak_comm"),
                    coef = coef(model.1.Mult.gof), statistics = c(rocpr), 
                    parallel = "multicore", ncpus = 3)
  
  gof.list.5[[i]] <- gof.ofs.E1
  
};note("Done.")

#saveRDS(gof.list.5, file="inference/gof/oos_model.5.rds")


















# --- Model 6 (Network effects + node effects for strong comms)
# NB: because bridging ties are not present in every year, we won't be able to estimate a coefficient during those absent years.
# so, either fit and estimate out-of-sample accuracy without the nodefactor("strong_bridge") term or subset the lists to
# only those years in which strong bridges are present. the results are substantively similar either way, so the models below 
# opt to retain the entire time range and skip the nodefactor("strong_bridge") term.

mid_sub_s <- readRDS("data/pauls_cranmer/prepped/subsets/strong_nodes/outcome_nets/mid_sub_s_weaknodes_prop0.2.rds")
directcont_sub_s <- readRDS("data/pauls_cranmer/prepped/subsets/strong_nodes/predictors/nodeset_weak/directcont_sub_s_weaknodes.rds")
capratio_sub_s <- readRDS("data/pauls_cranmer/prepped/subsets/strong_nodes/predictors/nodeset_weak/capratio_sub_s_weaknodes.rds")
tradedep_sub_s <- readRDS("data/pauls_cranmer/prepped/subsets/strong_nodes/predictors/nodeset_weak/tradedep_sub_s_weaknodes.rds")
igosec_sub_s <- readRDS("data/pauls_cranmer/prepped/subsets/strong_nodes/predictors/nodeset_weak/igosec_sub_s_weaknodes.rds")
igoecon_sub_s <- readRDS("data/pauls_cranmer/prepped/subsets/strong_nodes/predictors/nodeset_weak/igoecon_sub_s_weaknodes.rds")



gof.list.6 <- list()
for (i in 1:15){
  
  print(i)
  
  j=i+4
  k=i+5
  lst=i:j
  
  directcont.nw.gof1 <- directcont_sub_s[lst]
  capratio.gof1 <- capratio_sub_s[lst]
  tradedep.gof1 <- tradedep_sub_s[lst]
  igosec.gof1 <- igosec_sub_s[lst]
  igoecon.gof1 <- igoecon_sub_s[lst]
  
  set.seed(1912)
  model.1.Mult.gof <- btergm(
    mid_sub_s[lst] ~ 
      edges + 
      altkstar(2, fixed = TRUE) + 
      cycle(4) + 
      gwesp(0, fixed = TRUE) + 
      edgecov(directcont.nw.gof1) +
      edgecov(capratio.gof1) + 
      edgecov(tradedep.gof1) + 
      edgecov(igosec.gof1) + 
      edgecov(igoecon.gof1) +
      nodematch("strong_comm"), 
    R = 1000
    , parallel = "multicore", ncpus = 3
  )
  
  
  directcont.nw.gof1 <- directcont_sub_s[k]
  capratio.gof1 <- capratio_sub_s[k]
  tradedep.gof1 <- tradedep_sub_s[k]
  igosec.gof1 <- igosec_sub_s[k]
  igoecon.gof1 <- igoecon_sub_s[k]
  
  set.seed(1912)
  gof.ofs.E1 <- gof(model.1.Mult.gof, nsim=50, target = mid_sub_s[[k]],
                    formula = mid_sub_s[k] ~ 
                      edges + 
                      altkstar(2, fixed = TRUE) + 
                      cycle(4) + 
                      gwesp(0, fixed = TRUE) + 
                      edgecov(directcont.nw.gof1) +
                      edgecov(capratio.gof1) + 
                      edgecov(tradedep.gof1) + 
                      edgecov(igosec.gof1) + 
                      edgecov(igoecon.gof1) +
                      nodematch("strong_comm"), 
                    coef = coef(model.1.Mult.gof), statistics = c(rocpr), 
                    parallel = "multicore", ncpus = 3)
  
  gof.list.6[[i]] <- gof.ofs.E1
  
};note("Done.")

#saveRDS(gof.list.6, file="inference/gof/oos_model.6.rds")















# --- Check AUCs--- #

gof.list.1 <- readRDS("inference/gof/oos_model.1.rds")
gof.list.2 <- readRDS("inference/gof/oos_model.2.rds")
gof.list.3 <- readRDS("inference/gof/oos_model.3.rds")
gof.list.4 <- readRDS("inference/gof/oos_model.4.rds")
gof.list.5 <- readRDS("inference/gof/oos_model.5.rds")
gof.list.6 <- readRDS("inference/gof/oos_model.6.rds")


mod.1.pr <- data.frame()
for(i in 1:15){
  gof.i <- gof.list.1[[i]]
  pr.i <- gof.i$`Tie prediction`[[4]]
  pr.r <- gof.i$`Tie prediction`[[5]]
  df <- data.frame(pr.i, pr.r)
  mod.1.pr <- rbind(mod.1.pr, df)
}
colnames(mod.1.pr) <- c("Model 1: Baseline", "Model 1: Random")


mod.2.pr <- data.frame()
for(i in 1:15){
  gof.i <- gof.list.2[[i]]
  pr.i <- gof.i$`Tie prediction`[[4]]
  pr.r <- gof.i$`Tie prediction`[[5]]
  df <- data.frame(pr.i, pr.r)
  mod.2.pr <- rbind(mod.2.pr, df)
}
colnames(mod.2.pr) <- c("Model 2: Weak Comms", "Model 2: Random")


mod.3.pr <- data.frame()
for(i in 1:15){
  gof.i <- gof.list.3[[i]]
  pr.i <- gof.i$`Tie prediction`[[4]]
  pr.r <- gof.i$`Tie prediction`[[5]]
  df <- data.frame(pr.i, pr.r)
  mod.3.pr <- rbind(mod.3.pr, df)
}
colnames(mod.3.pr) <- c("Model 3: Strong Comms", "Model 3: Random")


mod.4.pr <- data.frame()
for(i in 1:15){
  gof.i <- gof.list.4[[i]]
  pr.i <- gof.i$`Tie prediction`[[4]]
  pr.r <- gof.i$`Tie prediction`[[5]]
  df <- data.frame(pr.i, pr.r)
  mod.4.pr <- rbind(mod.4.pr, df)
}
colnames(mod.4.pr) <- c("Model 4: Strong Comms no contig", "Model 4: Random")


mod.5.pr <- data.frame()
for(i in 1:15){
  gof.i <- gof.list.5[[i]]
  pr.i <- gof.i$`Tie prediction`[[4]]
  pr.r <- gof.i$`Tie prediction`[[5]]
  df <- data.frame(pr.i, pr.r)
  mod.5.pr <- rbind(mod.5.pr, df)
}
colnames(mod.5.pr) <- c("Model 5: Weak", "Model 5: Random")


mod.6.pr <- data.frame()
for(i in 1:length(gof.list.6)){
  gof.i <- gof.list.6[[i]]
  pr.i <- gof.i$`Tie prediction`[[4]]
  pr.r <- gof.i$`Tie prediction`[[5]]
  df <- data.frame(pr.i, pr.r)
  mod.6.pr <- rbind(mod.6.pr, df)
}
colnames(mod.6.pr) <- c("Model 6: Strong", "Model 6: Random")






# --- AUC boxplots --- #

randoms <- c(mod.1.pr$`Model 1: Random`, mod.2.pr$`Model 2: Random`,
             mod.3.pr$`Model 3: Random`, mod.4.pr$`Model 4: Random`,
             mod.4.pr$`Model 5: Random`, mod.6.pr$`Model 6: Random`)


dt <- data.table("Model 1:\nBaseline" = mod.1.pr$`Model 1: Baseline`, 
                 "Model 2:\nWeak Comms" = mod.2.pr$`Model 2: Weak Comms`, 
                 "Model 3:\nStrong Comms" = mod.3.pr$`Model 3: Strong Comms`, 
                 "Model 4:\nStrong Comms\n(No Contig.)" = mod.4.pr$`Model 4: Strong Comms no contig`, 
                 "Model 5:\nWeak Node\nEffects" = mod.5.pr$`Model 5: Weak`,
                 "Model 6:\nStrong Node\nEffects" = mod.6.pr$`Model 6: Strong`)

dt2 <- data.table("Random Graph" = randoms)
dt.m <- melt(dt)
dt.m2 <- merge(dt.m, melt(dt2), all = T)


ggplot(dt.m2, aes(variable, value)) + 
  geom_boxplot(fill="gray40", outlier.colour = "gray50") +
  stat_summary(geom = "crossbar", width=0.65, fatten=0, color="white", 
               fun.data = function(x){c(y=median(x), ymin=median(x), ymax=median(x))}) +
  theme_minimal() + 
  theme(axis.text=element_text(size=16), axis.title.y = element_text(size = 16)) +
  labs(x = "", y = "AUC (PR)")

#ggsave("box_pr_mult_f.pdf", width=13, height=5)

